*************************************************************
* The following code is part of the replication archive of  * 
* Genovese, Kern & Martin (2016) "Policy Alternation", ISQ  * 
*************************************************************

/* The code replicates additional analyses reported in the Appendix of Genovese et al (2016). 
 Note: The Stata MSTAR code below is discussed in Franzese, Hays and Kachi (2009). */

***********************************************
* All distances: last 3 columns in Table A.5  *
***********************************************

* MSTAR model with 3 spatial matrices: geographical distance, IGOs and  trade dyads + subsidy 

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1 rho2 rho3 sigma
tempvar A rSL1 rSL2 rSL3
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
gen double `rSL3'=`rho3'*SL3
scalar p1 = `rho1'
scalar p2 = `rho2'
scalar p3 = `rho3'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix p3W3 = p3*W3
matrix IpW = I_n - p1W1 - p2W2 - p3W3
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'- `rSL3'-`mu', 0, `sigma'))
end


**********************
*Open Data For Weights
**********************

spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag.dta", name(W1) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/trade_in2_lag.dta", name(W2) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/IGO_en_in2_lag.dta", name(W3) standardize 


**************************
*Open Data for Regression
**************************

drop _all


use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear

global nobs = 233 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y greentax_pc_cons
global X lag_gtax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp l_subsidy_rgdppc country_2-country_27  

mkmat $Y, matrix(Y)

matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W2*Y
svmat SL2, n(SL2)
matrix SL3 = W3*Y
svmat SL3, n(SL3)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X) (rho1:) (rho2:) (rho3:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo D1

drop SL11
drop SL21
drop SL31

**********************************************

* MSTAR model with spatially lagged subsidy

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1 rho2 rho3 sigma
tempvar A rSL1 rSL2 rSL3
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
gen double `rSL3'=`rho3'*SL3
scalar p1 = `rho1'
scalar p2 = `rho2'
scalar p3 = `rho3'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix p3W3 = p3*W3
matrix IpW = I_n - p1W1 - p2W2 - p3W3
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'- `rSL3'-`mu', 0, `sigma'))
end


**********************
*Open Data For Weights
**********************

spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag.dta", name(W1) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/trade_in2_lag.dta", name(W2) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/IGO_en_in2_lag.dta", name(W3) standardize 


**************************
*Open Data for Regression
**************************

drop _all

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear

global nobs = 233 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y greentax_pc_cons
global X lag_gtax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp l_subsidy_rgdppc country_2-country_27  
global Z l_subsidy_rgdppc

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)

matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W2*Y
svmat SL2, n(SL2)
matrix SL3 = W3*Y
svmat SL3, n(SL3)
matrix SL4 = W1*Z
svmat SL4, n(SL4)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X SL4) (rho1:) (rho2:) (rho3:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo D2

drop SL11
drop SL21
drop SL31

**********************************************

* MSTAR model with interaction with spatially lagged subsidyXeconomic integration

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1 rho2 rho3 sigma
tempvar A rSL1 rSL2 rSL3
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
gen double `rSL3'=`rho3'*SL3
scalar p1 = `rho1'
scalar p2 = `rho2'
scalar p3 = `rho3'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix p3W3 = p3*W3
matrix IpW = I_n - p1W1 - p2W2 - p3W3
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'- `rSL3'-`mu', 0, `sigma'))
end


**********************
*Open Data For Weights
**********************

spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag.dta", name(W1) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/trade_in2_lag.dta", name(W2) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/IGO_en_in2_lag.dta", name(W3) standardize 

**************************
*Open Data for Regression
**************************

drop _all

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear 

global nobs = 233 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y greentax_pc_cons
global X lag_gtax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp l_subsidy_rgdppc SL4 country_2-country_27  
global Z l_subsidy_rgdppc

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)


matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W2*Y
svmat SL2, n(SL2)
matrix SL3 = W3*Y
svmat SL3, n(SL3)
matrix SL4 = W1*Z
svmat SL4, n(SL4)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

*ml model lf splag_ll (mu: $Y=$X c.SL4##c.actual_economic_flows) (rho1:) (rho2:) (rho3:) (sigma:)
ml model lf splag_ll (mu: $Y=$X c.SL4#c.actual_economic_flows) (rho1:) (rho2:) (rho3:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo D3

drop SL11
drop SL21
drop SL31


**************
* Figure A.4 * 
**************

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1 rho2 rho3 sigma
tempvar A rSL1 rSL2 rSL3
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
gen double `rSL3'=`rho3'*SL3
scalar p1 = `rho1'
scalar p2 = `rho2'
scalar p3 = `rho3'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix p3W3 = p3*W3
matrix IpW = I_n - p1W1 - p2W2 - p3W3
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'- `rSL3'-`mu', 0, `sigma'))
end


**********************
*Open Data For Weights
**********************

spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag.dta", name(W1) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/trade_in2_lag.dta", name(W2) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/IGO_en_in2_lag.dta", name(W3) standardize 


**************************
*Open Data for Regression
**************************

drop _all

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear 

global nobs = 233 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y greentax_pc_cons
global X lag_gtax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp l_subsidy_rgdppc  country_2-country_27  
global Z l_subsidy_rgdppc

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)


matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W2*Y
svmat SL2, n(SL2)
matrix SL3 = W3*Y
svmat SL3, n(SL3)
matrix SL4 = W1*Z
svmat SL4, n(SL4)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X c.SL4##c.actual_economic_flows) (rho1:) (rho2:) (rho3:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

margins, at(SL4 =(0(25)8000) actual_economic_flows =(30)) post
estimates store marginplot3

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1 rho2 rho3 sigma
tempvar A rSL1 rSL2 rSL3
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
gen double `rSL3'=`rho3'*SL3
scalar p1 = `rho1'
scalar p2 = `rho2'
scalar p3 = `rho3'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix p3W3 = p3*W3
matrix IpW = I_n - p1W1 - p2W2 - p3W3
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'- `rSL3'-`mu', 0, `sigma'))
end

**********************
*Open Data For Weights
**********************

spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag.dta", name(W1) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/trade_in2_lag.dta", name(W2) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/IGO_en_in2_lag.dta", name(W3) standardize 

**************************
*Open Data for Regression
**************************

drop _all

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear 

global nobs = 233 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y greentax_pc_cons
global X lag_gtax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp l_subsidy_rgdppc  country_2-country_27  
global Z l_subsidy_rgdppc

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)


matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W2*Y
svmat SL2, n(SL2)
matrix SL3 = W3*Y
svmat SL3, n(SL3)
matrix SL4 = W1*Z
svmat SL4, n(SL4)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X c.SL4##c.actual_economic_flows) (rho1:) (rho2:) (rho3:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

margins, at(SL4 =(0(25)8000) actual_economic_flows =(95)) post
estimates store marginplot4

cd "/Users/.../gkm_replication_isq/green taxes and subsidies/Figures"

graph set window fontface "Times New Roman"

coefplot  (marginplot3,  lpattern(ldash) lcolor(gs13))  , at ytitle(Marginal Effect of Economic Flows) xtitle(Spatially Weighted Green Subsidies) title("Low Levels of Economic Flows") recast(line) lwidth(*2) ciopts(recast(rline) lpattern(dash)) ylabel(0(200)800,nogrid) xsca(alt) yline(0, lcolor(red)) plotregion(lcolor(black)) scheme(s2mono) graphregion(fcolor(white))   saving(marginplot3, replace)
coefplot  (marginplot4,  lpattern(ldash) lcolor(gs13))  , at ytitle(Marginal Effect of Economic Flows) xtitle(Spatially Weighted Green Subsidies) title("High Levels of Economic Flows") recast(line) lwidth(*2) ciopts(recast(rline) lpattern(dash)) ylabel(0(200)800,nogrid) xsca(alt) yline(0, lcolor(red)) plotregion(lcolor(black)) scheme(s2mono) graphregion(fcolor(white)) saving(marginplot4, replace) 

replace  SL4=SL4*2
twoway histogram  SL4 if actual_economic_flows<75, freq gap(5) color(grey14) ylabel(0(20)40,nogrid) plotregion(lcolor(black)) xlabel(0(2000)8000,nogrid) fysize(22)  xtitle(Spatially Weighted Green Subsidies)  graphregion(fcolor(white)) saving(hist3, replace) 
twoway histogram  SL4 if actual_economic_flows>75, freq gap(5) color(grey14)  ylabel(0(20)40,nogrid) plotregion(lcolor(black))  xlabel(0(2000)8000,nogrid) fysize(22) xtitle(Spatially Weighted Green Subsidies)  graphregion(fcolor(white))  saving(hist4, replace) 

graph combine marginplot3.gph hist3.gph, hole(2 4) imargin(0 0 1 0) graphregion(margin(l=10 r=10)) graphregion(fcolor(white))  saving(pic3, replace)
graph combine marginplot4.gph hist4.gph, hole(2 4) imargin(0 0 1 0) graphregion(margin(l=10 r=10)) graphregion(fcolor(white))  saving(pic4, replace)  

gr combine pic3.gph pic4.gph, ycommon col(2) plotregion(fcolor(white)) graphregion(fcolor(white))   title(`"Dependent Variable: Green Tax Levels"', color(black)) saving(marginplot2and3, replace)

